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I ' Context. BL Lac objects undergo strong flux variations involving considerable changes in their spectral shapes. We specifically investigate the 
' X-ray spectral evolution of Mrk 421 over a time span of about nine years. 
^-H Aims. We aim at statistically describing and physically understanding the large spectral changes in X rays observed in Mrk 421 over this time 
span. 

Metiiods. We perform a homogeneous spectral analysis of a wide data set including archived observations with ASCA, BeppoSAX, RXTE, as 
well as published and unpublished XMM-Newton data. The presence of uncertainties is taken into account in our correlation analysis. The 
significance of the correlations found and possible spurious effects are studied with Monte Carlo simulations. 

Resuits. We find that the Mrk 421 spectral energy distribution (SED) has a lower peak at energies that vary in the range, 0.1-10 keV while its 
O I X-ray spectrum is definitely curved. Parameterizing the X-ray spectra with a log-parabolic model, we find a positive correlation between the 
position and the height of the SED peak. In addition, we find a negative trend of the spectral curvature parameter vs. the SED peak energy. 
Conciusions. We show that these relations between the spectral parameters are consistent with statistical or stochastic acceleration of the 
emitting particles, and provide insight into the physical processes occurring in BL Lac nuclei. 
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H 1. Introduction bump in the IR-to-optical band, as opposed to the HBL (High 

■ - - energy peaked BL Lac) that peak in the UV-X-ray band. 

In the popular unification scenario of the Active Galactic According to the Synchrotron Self Compton (SSC) emission 

Nuclei (AGNs) BL Lac objects are understood as sources in mechanism, the high frequency bump is to be attributed to in- 

which a relativistic jet is produced by the central engine and ^^^^^ Compton scattering of synchrotron photons by the same 

points close to our hne of sight, see Antonucci (1993) and Urry p^^^i^^^^^ of relativistic electrons that produce the synchrotron 

& Padovam (1995). These objects are marked by featureless ej^ission (Jones, O'Dell & Stein 1974, GhiselUni & Maraschi 

spectra and a high flux from the radio to the y-ray band, en- ^939). Both emissions are produced in a relativistic outflow 

dowed with strong variability and with high polarisation ob- ^^^j^ ^ ^^j^ Lorentz factor F ^ 10; when observed at an an- 

served in the optical band. In fact, the emission from BL Lacs ^ ^j^^y subjected to the eff^ects of a beaming factor 

constitutes one of the best examples of direct non-thermal ra- g - i /(r(l -yScos 0)) 
diation from particles not in thermal equilibrium with photons 

nor among themselves. So these sources are keenly interesting With its redshift z = 0.031, Mrk 421 is among the closest 

to investigate particles acceleration mechanisms in AGNs. and best studied HBL. In fact, it is one of brightest BL Lac 

On a close look, their Spectral Energy Distribution (SED) objects in the UV and in the X-ray bands, observed in y rays 

appears as double-peaked. The lower energy component is EGRET (Lin et al. 1992); it was also the first extragalactic 

widely interpreted as synchrotron emission from highly rela- source detected at TeV energies in the range 0.5-1.5 TeV by the 

tivistic electi-ons with Lorentz factors y in excess of 10^; of- Whipple telescopes (Punch et al. 1992, Petry et al. 1996). 

ten it peaks at frequencies in the IR to the X-ray band. The The source is classified as HBL because its synchrotron 

actual position of this peak has been suggested by Padovani emission peak ranges from a fraction of a keV to several keVs. 

& Giommi (1995) as a marker for a classification; they de- In fact, its flux changes go along with strong spectral varia- 

fine LBL (Low energy peaked BL Lac) objects with the first tions (Fossati et al. 2000a). The spectral shape generaUy ex- 

hibits a marked curvature, well described by a log-parabolic 

Send ojfprint requests to: andrea.tramacere@romal.infn.it model (Landau et al. 1986). 
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We present here a new study of the time- variable properties 
of its SED based on the full collection of XMM-Newton data 
taken with EPIC CCD cameras, sensitive in the energy range 
0.3-10 keV; we also use previous X-ray observations with 
ASCA, BeppoSAX, RXTE, EUVE to cover an overall time 
span of nine years. In particular, XMM-Newton observed 
Mrk 421 many times, and a large amount of data has been col- 
lected. Some of these observations have been recently analysed 
by other authors, but a full report of all observations between 
May 2000 and November 2005 had not been published yet. 
We perform a new analysis of all XMM-Newton observations 
of Mrk 421, particularly relevant as its synchrotron emission 
peaks in the XMM-Newton energy range. 

Specifically, extending the previous work of Tanihata et 
al. (2004), we study correlations between the position and the 
height of SED peak, and interpret them in terms of signatures 
of the synchrotron emission. We also investigate a new corre- 
lation of the SED peak with the spectral curvature, and discuss 
its impUcations as to the electron acceleration mechanisms. 

This paper is organized as follows. Sect. 2 describes the 
data set and our procedure to reduce the XMM-Newton data. 
Sect. 3 reports our data analysis, including temporal and spec- 
tral investigations. The statistical analysis and the comparison 
with numerical simulations and theoretical interpretations are 
discussed in Sect. 4 and Sect. 5. Our conclusions are given and 
discussed in Sect. 6. 

2. Data sets and Reductions 

2.1. XMM-Newton observations 

Mrk 421 was observed by XMM-Newton between 25/05/00 
and 9/11/05, on more than twenty occasions, by means of 
all EPIC CCD cameras: the EPIC-PN (Struder et al. 2001), 
and EPIC-MOS (Turner et al. 2001), operating in different 
modes and with different filters as described in Table 1. 
These data were reduced as follows. Extractions of all Ught 
curves, source and background spectra were done using the 
XMM-Newton Science Analysis System (SAS) v6.5.0. The 
CaUbration Index File (CIF) and the summary file of the 
Observation Data File (ODF) were generated using Updated 
Calibration File (CCF) following the "User's Guide to the 
XMM-Newton Science Analysis Syste" (Issue 3.1) (Loiseau 
et al. 2004) and "The XMM-Newton ABC Guide" (vers. 2.01) 
(Snowden et al. 2004). Event files were producted by XMM- 
Newton EMCHAIN pipeUne. The standard reduction of the 
events list for MOS data, was performed involving subtrac- 
tion of hot and dead pixels, removal of events due to the elec- 
tronic noise and correction of event energies for charge transfer 
losses. 

To provide the most conservative screening criteria MOS 
data files were also filtered to include all single to quadruple 
events (PATTERN < 12) with pulse high rate in the range of 
500 to 12,000 eV and with expression FLAG = 0. Lightcurves 
for every dataset were extracted and all high-background time 
intervals, were filtered out by excluding time interval contami- 
nated by solar flare signal. To perform this selection, the count 
rate in the 10 - 15 keV band for the entire MOS detectors were 



Table 1. XMM-Newton Log of observations of Mrk 421 . 



Date 




EPIC-MOS 1 




EPIC-MOS2 


dd/nun/yy 


Frame 


rlllcr 


Exp 


Frame 


rlllcr 


Exp 


25/05/2000 


FU 


Md 


18000 


PW 


Md 


24001 


01/05/2000 


PW 


Md 


37251 


PW 


Md 


37250 


13/05/2000 


PW 


Md 


46950 


PW 


Md 


46951 


14/05/2000 


PW 


Md 


41948 


PW 


Md 


41948 


08/05/2001 


PW 


Tn 


38408 


PW 


Tn 


38409 


04/05/2002 


FU 


Tn 


38898 


FW 


Tn 


39163 


05/05/2002 


FW 


Tk 


19422 


FW 


Tn 


19419 


04/05/2002 


FU 


Tn 


23407 


FU 


Tn 


23409 


04/05/2002 


FU 


Tn 


23416 


FU 


Tn 


23410 


14/05/2002 


FW 


Tn 


71200 


FW 


Tn 


71200 


14/05/2002 


FW 


Tn 


11172 


FW 


Tn 


11174 


15/05/2002 


FW 


Tn 


11175 


FW 


Tn 


11174 


01/05/2002 


FW 


Tn 


70873 


FW 


Tn 


70873 


02/05/2002 


FW 


Tn 


11173 


FW 


Tn 


11172 


02/05/2002 


FW 


Tn 


11171 


FW 


Tn 


11175 


01/05/2003 


FU 


Md 


40542 


PW 


Md 


40923 


02/05/2003 


FU 


Tk 


8912 


PW 


Tk 


8911 


07/05/2003 


FU 


Tk 


8000 


PW 


Tk 


8000 


14/05/2003 


PW 


Tn 


48668 


FU 


Tn 


48413 


10/05/2003 


PW 


Md 


25848 


PW 


Tk 


25851 


06/05/2004 


PW 


Md 


61532 


PW 


Md 


61534 


07/05/2005 


RF 


Tn 


9765 


RF 


Tn 


9771 


07/05/2005 


RF 


Md 


9465 


RF 


Md 


9470 


09/05/2005 


PW 


Tk 


59809 


PW 


Tk 


59802 



(a) FRAME: PW=partial window, FU=fast uncompressed, FW=full 
window, RF=refresh frame. 

(b) FILTER: Tn=thin, Md=Medium, Tk=thick. 



determined. We first discarded only time intervals with a count 
rate that exceeds 0.35 counts per second as indicated in the 
"User's Guide to the XMM-Newton Science Analysis System" 
(Issue 3.1) (Loiseau et al. 2004). However applying this crite- 
rion, we noticed that a low background state, placed between 
two high neighboring peaks due to a solar flare, can even in- 
clude a residual contamination to the source signal which mod- 
ifies the spectral distribution. We adopted a more conserva- 
tive selection for good time intervals, because we preferred 
to have a high quality signal, excluding time ranges that ap- 
peared contaminated by solar flares. Then we selected good 
time intervals by direct inspection far from solar flare peaks 
and without count rate variations on time scales shorter than 
500 seconds. The TABGTIGEN task of XMM-Newton Science 
Analysis System (SAS) was used to build good time intervals. 

Photons were extracted from an annular region using dif- 
ferent apertures to minimize pile-up, which affects MOS data. 
The mean value of external radius for the annular region is 40". 
To filter out pixels affected by significant pile-up, the internal 
region was selected by using EPATPLOT task in XMM-Newton 
Science Analysis System (SAS) for each observation. 

In FULL WINDOW images, the background spectrum 
was extracted from a circular region of the size compara- 
ble to the source region, in a place where visible sources 
were not present (typically off axis). For other observations, 
in PARTIAL WINDOW images, no regions sufficiently far 
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from the source for the background extraction were found. In 
these cases we used background from blank-field event files 
(www.sr.bham.ac.uk). Anyway we estimated that the average 
X-rays background flux was always at ~ 1% level of source 
flux, resulting in a negligible contamination in the spectral pa- 
rameter determination. The Photon Redistribution Matrix and 
the Ancillary Region File were created for each observation, 
by using RMFGEN and ARFGEN tasks of SAS. A more re- 
stricted energy range (0.5-10 keV) was used to account for 
possible residual caUbration uncertainties. To insure the valid- 
ity of Gaussian statistics, data were grouped by combining in- 
strumental channels so that each new bin comprises 40 counts 
or more. 

2.2. Other X-ray observations 

To study the spectral behaviour of Mrk 42 1 over a wide time 
interval we also used other X-ray observations performed with 
ASCA, RXTE, EUVE and BeppoSAX, with references given 
below. All these observations were analysed in terms of the 
same spectral model as the XMM-Newton ones homogeneous 
(described in §3), so we have obtained a sample of spectral 
parameters that covers about 9 years of Mrk 421 observations 
in X rays. 

Mrk 421 was pointed continuously by ASCA for 7 days 
between 23 April of 1998 and 30 April 1998, with all its detec- 
tors on: two SISs (Stage Imaging Spectrometers) and two GISs 
(Gas Imaging Spectrometers). To extend tiie ASCA data base to 
higher energy, these observations were coordinated with RXTE, 
specifically with the PCA (Proportional Counter Array; 2-60 
keV) and the HEXTE (High Energy X-ray Timing Experiment; 
15-200 keV) instruments on board RXTE. We refer to Tanihata 
et al. (2001) for the ASCA data reduction and spectral anal- 
ysis, and to Tanihata et al. (2004) for the the combination of 
RXTE with A5CA. 

BeppoSAX observed Mrk 421 on 7 occasions dur- 
ing 1997, in which all instruments on board were operat- 
ing: LEGS (Low Energy Concentrator Spectrometer; 0.1-10 
keV), MECS (Medium Energy Concentrator Spectrometer; 
1.3-10 keV) and the PDS (Phoswich Detector System; 13-300 
keV). Spectral analyses of these observations are reported in 
Massaro et al. (2004) and references therein. Concerning 1998 
BeppoSAX observations of Mrk 421, started 3 days before the 
ASCA pointing, Tanihata et al. (2004) reported the spectral 
analysis of data obtained by the LEGS and the MECS (see also: 
Maraschi et al. 1999, Fossati et al. 2000a, 2000b, Massaro et al. 
2004) . In May 1 999 Mrk 42 1 was pointed for about 4 day s ; two 
long uninterrupted observations were performed from 26 April 
and 3 May, 2000, followed by another observation from 9 May 
2000 and 12 May 2000. In all observations of 1999 and 2000, 
Mrk 421 showed a strong variability on time scales of days. 
Data reduction and spectral analysis of these BeppoSAX ob- 
servations are described in the references above. 

3. Spectral Analysis 

Our XMM-Newton spectral analysis was performed with the 
xsPEc software package, version 11.3 (Amaud, 1996). Before 



Table 2. Spectral parameters of the LP model of time resolve 
spectra for the 0.5-10 XMM-Newton EPIC-MOS observations 
of Mrk 421 



Date E], b 5;/ xUd-O-f-) 



01/11/00 


0.22 (0.08) 


0.25 


(0.05) 


132.8 (9.6) 


0.86 


(174) 




0.22 (0.04) 


0.23 


(0.02) 


124.8 (4.8) 


0.95 


(307) 




0.31 (0.06) 


0.28 


(0.03) 


136.0 (4.8) 


1.12 


(225) 


13/11/00 


1.09 (0.04) 


0.27 


(0.02) 


340.6 (1.0) 


1.06 


(372) 




0.92 (0.03) 


0.40 


(0.02) 


342.9 (1.1) 


1.07 


(363) 


08/05/01 


0.74 (0.05) 


0.35 


(0.02) 


254.4 (1.6) 


1.01 


(285) 




0.82 (0.04) 


0.27 


(0.02) 


263.5 (11) 


1.23 


(347) 




0.85 (0.04) 


0.37 


(0.02) 


249.4 (1.1) 


1.11 


(311) 


04/05/02 


0.39 (0.07) 


0.49 


(0.06) 


147.2 (6.4) 


0.98 


(139) 




0.38 (0.06) 


0.46 


(0.05) 


144.0 (6.4) 


1.06 


(154) 




0.44 (0.06) 


0.56 


(0.06) 


120.0 (4.8) 


0.96 


(153) 


05/05/02 


0.23 (0.04) 


0.44 


(0.04) 


102.4 (6.4) 


0.98 


(186) 


14/11/02 


1.43 (0.04) 


0.39 


(0.03) 


347.2 (1.6) 


1.00 


(221) 




2.0 (0.1) 


0.32 


(0.05) 


380.8 (4.8) 


1.22 


(155) 


01/12/02 


0.87 (0.05) 


0.41 


(0.03) 


172.6 (1.1) 


1.01 


(228) 




0.72 (0.07) 


0.38 


(0.04) 


169.6 (1.2) 


1. 10 


(195) 




0.67 (0.06) 


0.40 


(0.04) 


167.7 (1.4) 


0.91 


(212) 


02/12/02 


0.64 (0.05) 


0.34 


(0.03) 


196.8 (1.6) 


1.00 


(216) 




0.83 (0.08) 


0.35 


(0.04) 


227.2 (1.6) 


1.08 


(190) 


02/06/03 


0.36 (0.08) 


0.45 


(0.06) 


51.2 (3.2) 


1.09 


(139) 


14/11/03 


1.20 (0.03) 


0.43 


(0.02) 


467.2 (1.6) 


1.15 


(321) 




0.79 (0.04) 


0.48 


(0.03) 


441.6 (3.2) 


1.10 


(240) 


10/12/03 


0.87 (0.05) 


0.38 


(0.03) 


209.1 (1.1) 


1.04 


(274) 




0.89 (0.05) 


0.42 


(0.03) 


207.2 (1.3) 


1.17 


(240) 


06/05/04 


1.74 (0.04) 


0.33 


(0.03) 


432.0 (1.6) 


1.09 


(332) 




2.4 (0.1) 


0.24 


(0.02) 


494.4 (3.2) 


1.24 


(335) 


07/11/05 


0.9 (0.1) 


0.38 


(0.06) 


284.8 (3.2) 


1.06 


(145) 




0.8(0.1) 


0.44 


(0.07) 


296.0 (4.8) 


1.04 


(125) 


09/11/05 


0.63 (0.03) 


0.40 


(0.02) 


420.8 (3.2) 


1.27 


(365) 




0.71 (0.04) 


0.41 


(0.02) 


430.4 (3.2) 


1.25 


(302) 




0.68 (0.04) 


0.36 


(0.02) 


460.8 (3.2) 


1.16 


(300) 




0.79 (0.05) 


0.40 


(0.03) 


505.6 (3.2) 


1.22 


(266) 



(*) Ep in keV 

(**)Spin IQ-'^ergcm-^ s"' 



evaluating any spectral curvature in the 0.5-10 keV range, it 
is important to reliably assess the absorption at low energies 
due to interstellar gas. In the optical, high resolution images 
of the host, early-type galaxy of Mrk 421 do not show in the 
brightness profile any evidence of large amounts of absorbing 
material (Urry et al. 2000). In X-rays, it has been shown by 
Fossati et al. (2000b) and confirmed by Tanihata et al. (2004) 
and Massaro et al. (2004) that describing the spectral shape in 
terms of absorption not only would require a column density 
much higher than the Galactic value Nh = 1.61 x lO^^cm"^ 
(Lockman & Savage, 1995), but also would yield in any case 
unacceptable fits with very high;^'^- These findings motivate us 
to consider intrinsicaUy curved spectra with the column den- 
sity fixed at the Galactic value, in agreement with the above 
analyses of Mrk 421. 

In fact, we tried a number of different spectral models: 
single (PL) and broken (BPL) power-laws; a power-law with 
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Fig. 1. XMM-Newton EPIC-MOS spectrum from the observation of Mkn 421 performed on 08/05/01. Left: the systematic 
deviations on both sides of the residuals from a best fit PL with Galactic A^h show the need of intrinsic curvature. Right: the 
deviations disappear with the LP model with Galactic Nu ■ 



an exponential cutoff at high energy (PLC); and a logarithmic 
parabola (LP) in the form 

FiE) = K E-"-'' E ^ 

where a is the spectral slope at £1 - 1 keV, and b measures the 
curvature. 

We find that the PL model is not adequate to describe the 
spectral shape, consistent with other X-ray analyses reported 
above. In our observations we mostly obtained unacceptable 
values forxl (higher than 1.5); in the few cases where the 
is acceptable, the residuals show systematic deviations clearly 
indicating the presence of a curvature, see Fig. 1. PLC and 
BPL models give better spectral descriptions, but they also have 
drawbacks. The PLC model often estimates the cut-off energy 
at much higher than 10 keV, indicating that in the instrumental 
range there is a mild curvature rather than a sharp exponential 
cutoff. The BPL model gives systematic excess at high ener- 
gies. The LP model describes the spectral shape more satisfac- 
torily than all the above ones, giving systematically lower 
values and fewer residuals. 

The main goal of our analysis is the search of possible cor- 
relations between the spectral curvature (b), the SED peak en- 
ergy (Ep), and the corresponding SED peak value (S p). To esti- 
mate the correlation between the parameters of the spectral en- 
ergy distribution S(E) = E^F(E) we often use in place of Eq. 
(1) the equivalent functional form (LPS), see Sect. 3 in Massaro 
et al. 2004 and Tanihata et al. 2004 

5(£) = 5p lO-^'os^*^/^"', (2) 

where S p = Ep F(Ep). In this form the values of the parame- 
ters b, Ep and 5 p are estimated in the fitting procedure, whereas 
those are derived from Eq. (1) are affected by intrinsic correla- 
tions. For our analyses, we added a specific spectral modelling 
routine to the standard xspec v. 11.3 software package. 



We shall see that the estimates of the curvature parameter 
b are sensitive to changes of the source state during a pointing; 
so a time averaged spectrum may have a different curvature 
from the the time resolved ones. To avoid any such bias, we 
segmented our data at times when significant (3 sigma) count 
rate variations occur, subject to the constraint of having in each 
time resolved spectrum more than 120 spectral bins after rebin- 
ning as indicated in §2. The duration of our time segments was 
shorter than a few 10* s. 

An example of this effect is apparent from the two light 
curves plotted in Fig. 2 in the energy range 0.5-10 keV (left 
panel) and 4-10 keV (right panel) of the same observation 
on 13/11/00. The significant decrease in the count rate of the 
higher-energy curve is only barely detecable in the total light 
curve. Then the spectral analysis of the two time regions (A and 
B, see Fig. 2) shows a spectral dynamics not apparent when the 
whole data set is considered. 

The spectral best fit values for all XMM-Newton observa- 
tions are reported in Table 2. All statistical errors refer to the 
68% confidence level (one Gaussian standard deviation). 

Our results generally agree with those from observations 
with BeppoSAX, ASCA and Swift (see the references above 
and Tramacere et al. 2006 for the Swift analysis). In particu- 
lar, we find values of SED peak energy and of the curvature 
parameter consistent with their values, but we find differences 
with other XMM-Newton analyses. Specifically, Sembay et 
al. (2002) analysed XMM-Newton observations during 2000, 
2001 and 2002, but they used an energy range restricted to 2- 
10 keV and did not search for any curved shape. Ravasio et 
al. (2004), who studied EPIC-PN observations performed on 
04/1 1/02, 14/1 1/02 and 01/12/02, didn't found a significant im- 
provement using the log-paraboUc model respect to the power- 
law one. The difference between their and our analysis is likely 
due to our more conservative criterion to subtract solar flares in 
the data reduction. 
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Fig. 2. {Left panel): XMM-Newton light curve in the energy range 0.5-10 keV. (Right panel)'. XMM-Newton light curve in the 
energy range 4-10 keV, where count rate variations are evident between the two selected time regions. 



4. Testing the dynamics of synchrotron emission 

Correlations between S p and Ep provide interesting indica- 
tions concerning the driver of the spectral changes in X-rays, in 
terms of the synchrotron emission mechanism from one dom- 
inant homogeneous component. In this framework (Rybicki & 
Lightman, 1979) the dependence of 5 p on Ep is in the form of 
a power-law: 

Sp^E-p. (3) 

In fact, the synchrotron SED is expected to scale as 5 oc 
N 6^ at the energies E oc y^B 6, in terms of total emitter 
number A^, the magnetic field B, the typical electron energy 
ymc^, and the beaming factor 6, Thus a = 1.5 applies when the 
spectral changes are dominated by variations of the electron 
average energy; a = 2 as for changes of the magnetic field; 
a = 4 if changes in the beaming factor dominate; formally, 
or = 00, applies for changes only in the number of emitting 
particles. 

On a heuristic stand, inspection of Fig. 3 indicates that the 
two last cases are unlikely; next we proceed to a formal statis- 
tical analysis. 

4. 1. Correlation between Sp and Ep 

We investigate whether and how Sp and Ep correlate. 
The analysis is performed on two data sets: the XMM set 
and the FULL set . The former corresponds to the results 
of the spectral analysis in §3, while the latter is obtained on 
adding the results from previous analyses as reported in Sect. 3. 

We test for correlations between the spectral parameters 
on using both a linear correlation coefficient (run) and a log- 
arithmic one (riog), and fit the data with a Unear and log- 



linear model following the statistical methods described in 
D'Agostini (2005) (see also Appendix A). 

Results from our correlation study are given in Table 3. 
The correlation coefficient riog - 0.7 is significant, specifically 
with a chance probability value close to zero. The data in the 
XMM set may be described by a power-law function as in Eq. 
(3); the value of a so obtained is a = 1.4 + 0.3 . On fitting the 
FULL set we find the power-law index a = 1 .2 ± 0. 1 , consistent 
with the above value. Even if we took the range of Ep between 
0.5 - 1.0 keV where the run of the FULL data apparently 
steepens, the values of the best fit would yield a = 1.2 + 0.1, 
quite far from 4. This analysis confirms that the cases a = 1.5 
and or = 2 are those most relevant as dominant mechanisms. 

A closer analysis and an interpretation of these results are 
given in the next subsection. 

4.2. Signatures of tiie emission process 

A sensitive point in our statistical analysis and its interpretation 
is the presence of "hidden" parameters contributing to variabil- 
ity and/or to correlations. On the basis of the formahsm that for 
reader's convenience we recall in Appendix A, in the following 
we discuss in details the possibility that a correlation may be 
introduced by a parameter that is not directly determined from 
the data (whence the name "hidden"). 

This issue arises when one estimates two parameters, say x 
and y, statistically independent but depending physically on a 
third parameter, k say, that is not directly observed ("hidden"); 
in other words, we observe x = x(x', k) and y = y(y', k). Again, 
a correlation will arise between x and y due to their dependence 
on k, and this may add to any physical correlation between x' 
and/. 
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Fig. 3. The scatter plot of the peak energy Ep and the maximum of the SEDs S p compared with the results of a Monte Carlo sim- 
ulation (grey squares). Dashed line represents the power-law best fit without taking into account extravariance for the FULL data- 
set. Solid line represents the power-law best fit taking into account extravariance for the FULL data-set. 



Table 3. Statistical parameters for the correlation between Ep and S p. 

Dataset ru„ pn,, n„g pi„g a\. 

XMM 0.63 < 0.001 0.71 < 0.001 1.4 ±0.3 0.7 ±0.1 0.35 ± 0.05 

FULL 0.56 < 0.001 0.67 < 0.001 1.2 ±0.1 0.56 ± 0.05 0.33 ± 0.02 

FULL{Ep< 1.2) 0.72 < 0.001 0.72 < 0.001 1.7 ±0.1 0.99 ± 0.07 0.29 ± 0.02 

FULL(Ep> 1.2) 0.87 0.004 0.84 0.007 0.6 ±0.1 0.43 ± 0.08 0.13 ±0.03 



A case in point is provided by the beaming effects. In BL 
Lacs beaming is a key property; though only rarely measured 
directly, this can introduce dependence between the observed 
quantities S p and Ep. It is convenient to distinguish the vari- 
ables E'p, S'p in the beam reference frame from the observed 
ones expressed as 

Ep = e;s (4) 
Sp ^ S'p6'^. 

While this corresponds to the case a - 4- that our preliminary 
analysis tended to exclude, it is important to reconsider 6 
as a hidden parameter introducing a covariance term and 
contributing to correlation. 

Clearly, the actual contribution depends on the probability 
density functions (PDF) of Ep, S p and 6. To study this, we have 
generated with a Monte Carlo code a sample of uncorrected 
pairs of variables (x',y'). We then have transformed them ac- 



cording to Eq. (4), that is, with x = 6x' and y = 6 y', choosing 
6 = 10 as a typical value for Mrk 421. The x' variables were 
generated so that the x variables have the same PDF as Ep; on 
the other hand, the y' variables were generated from a normal 
distribution with mean value fi' - 1.0 and cr' - /j.' /3. We im- 
pose that the y variables have the same dispersion as 5 ^ in the 
FULL data; this dispersion depends not only on cr' but also on 
the distribution of S, so we have to derive a constraint on the 
value of (Tg. In fact, the standard deviation of the observed val- 
ues of S p is about 120 (erg cm"^ s"') to be compared with 
the average value of about 290 (Fig. 3); assuming, as an ex- 
treme case, that the dispersion cr^ of S p is generated only by 
the variance erg, and applying standard propagation, we obtain 
crs/fis^crs/(4Sp)< 10%. 

We now generate 6 from a normal distribution with fig - 
10 and cTg = 0.75, and find a coiTelation coefficient riog = 0.3 
between the logarithms of x and y. On decreasing the variance 
of y' and increasing that of 6, the correlation increases but only 
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slightly; for example, if we take cr' = /j'/6 and increase the 
value of (Tg to 0.95, we obtain ri„g - 0.36. 

In summary, the beaming factor can affect the observed 
correlation coefficient, but only up to values riog < 0.3 which 
are significantly below the observed value 0.67. We conclude 
that beaming alone is responsible neither for the values of 
a ^ I nor for the correlation observed. We stress that the value 
cr^/S < 0.1, bounded by the dispersion of the full data, is re- 
markably low when compared to the average typical value of 
/xg ^ 10, and imphes the beaming factor to have been closely 
constant for Mrk 421 during our observational span of about 9 
yr. 

It follows that both the variations of 5^ corresponding to 
the estimated values of a and their scatter must be importantly 
contributed by a physical process in the beam rest frame, such 
as variations of magnetic fields or scahng up or down of all the 
electron energies. 

In fact, a second eff'ect coming from a "hidden" parame- 
ter with its fluctuations is to produce scatter of the data. As 
explained in Appendix A, it is possible to account for this ef- 
fect by adding an "extravariance" in the hkehhood function. 
In the case of the XMM data set, the extravariance accounting 
for their scatter is estimated at cr^, = 0.35 + 0.05, and for the 
corresponding value of the slope we obtain ay = 0.7 + 0.1, 
whilst for the FULL data set we obtain cTv = 0.33 ± 0.02 and 
a^, = 0.56 + 0.05 (see Table 3). The derived values of the a in- 
dex are considerably affected by the extravariance; in fact, the 
extravariance term in the log-likelihood function dominates the 
term from measurement uncertainties, to the point of providing 
power-law indices close to those obtained simply from fitting 
the data with no weights for their precision. 

Such diff'erent values of a and ay (Tab. 3) lead us to test an 
actual change in the power-law index in Eq. (3) at about 1 keV. 
To do this, we split the FULL set into two sets, one with Ep < 
1 .2 keV and the other with Ep > 1.2 keV; results of fits are 
reported in the last two lines of Table 3. The values of a found 
for the two data sets diff'er significantly, being ori = 1.7±0.1and 
o'2 = 0.6+0.1, for Zip < 1.2keVand£'p > 1 .2 keV, respectively. 
This could be interpreted as a change of a with the state of the 
source, and we estimate the associated hkehhood by means of 
a Monte Carlo simulation. We generate a set of values for Ep, 
S p as follows: 

- Ep values are generated as to have the same PDF as ob- 
served once transformed by the first of Eq. (4). 

- Beaming factors are generated from a normal distribution 
with /^^ = 10 and as = 0.5 . 

- The slope of the (unbeamed) power-law is generated from 
two normal distributions with mean value 0.85 and standard 
deviation 0.09 for Ep < 1.2 keV, and mean value 0.33 and 
standard deviation 0.07 for Ep> 1.2 keV. 

- Values oi S p are generated from Eq. (3) and transformed 
by the second of Eq. (4). 

The number of events so generated is 150 (similar to the 
FULL data set), and in Fig. 2 we report their scatter plot (grey 
box). Interestingly, the simulated and the observed points have 
similar behaviour and similar statistical properties. We see that 
the correlations, given in the beam reference frame by ai and 



02, when subject to beaming with a narrow variance come very 
close to account for the data with their scatter. 



5. Probing signatures of the electron acceleration 

Finally, we perform a study of correlation between b and Ep, 
aimed at pinpointing possible signatures of the electron accel- 
eration processes in the spectral evolution of a homogeneous 
dominant component. The same data set and statistical tools of 
the previous section are used. 

5.1. Test of correlation between b and Ep 

Such a correlation had never been tested previously. An analy- 
sis of XMM set shows that the correlation coefficients are neg- 
ative and low, that is, r;,>, = -0.29 and riog - -0.13 (defined in 
Sect. 4). The former corresponds to a low chance occurrence 
(0.05), but this is significantly high (0.24) for riog. By direct 
inspection we see that there are three points (enclosed by the 
dashed ellipse in Fig. 4) corresponding to the same XMM point- 
ing on 01/1 1/00, that maximally deviate from the average sam- 
ple behaviour. On excising this observation and reanalysing the 
remaining data (denoted by XMM * in Table 5) we obtain the 
substantiaUy higher correlation coeflficients r;,„ = -0.60 and 
riog - -0.62, corresponding to the high significance, given in 
Table 4. 

To see whether the 01/11/00 observation constitutes a 
unique event, we analysed the FULL data set with and with- 
out (denoted by FULL *) this pointing; Table 4 shows that the 
values of r;,„ and r^g for FULL and FULL * are close, indicat- 
ing that the observation excised has not materially changed the 
overall statistical behaviour. 

These values are sufficiently high to warrant a brief dis- 
cussion of the importance of the correlation found, along two 
ways. This pointing may constitute just a rare source state in 
our sample. On the other hand, this may not constitute a sin- 
gular event, considering that all XMM and FULL data (Fig. 4) 
appear to outline an upper limit to the curvature observed at 
any values of Ep. For example, for Ep around 1 keV, we never 
observe curvatures higher than about b « 0.45. 

5.2. Interpretations 

The inverse correlation between b and Ep may be interpreted in 
the framework of acceleration processes of the emitting elec- 
trons. 

A first interpretation of the correlation between b and Ep is in 
the framework of statistical acceleration (Massaro et al. 2006 
and references therein). In this scenario the probability for a 
particle to be further accelerated decreases at high energies, be- 
ing inversely proportional to the energy itself. For example, this 
may occur when the particles are confined by a magnetic field, 
and the confinement efficiency decreases as the gyration radius 
increases. 
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Fig. 4. The scatter plot of the peak energy Ep and the curvature parameter b. The dashed line represents the power-law best fit 
without taking into account extravariance for the FULL data-set. The solid line represents the power-law best fit on taking into 
account extravariance for the FULL data-set. The dashed-dotted line represents the relation between Ep and b expected from Eq. 
(11) 



Table 4. Statistical parameters for the correlation between Ep and b 



Dataset 




Plin 




Phg 






cr,. 


XMM 


-0.29 


0.05 


-0.13 


0.24 


-0.04 ± 0.03 


-0.04 ± 0.06 


0.18 ±0.03 


XMM* 


-0.60 


< 0.001 


-0.60 


< 0.001 


-0.20 ± 0.03 


-0.22 ± 0.06 


0.12 ±0.02 


FULL 


-0.63 


< 0.001 


-0.70 


< 0.001 


-0.404 + 0.005 


-0.31 ±0.02 


0.13 ±0.01 


FULL* 


-0.67 


< 0.001 


-0.79 


< 0.001 


-0.404 ± 0.005 


-0.34 ± 0.02 


0.12 ±0.02 



In such cases the electron energy distribution is curved into 
a log-parabolic shape, and its curvature r is related to the frac- 
tional acceleration gain e as given by Massaro et. al. (2006): 

1 



log e 



(5) 



Note that r decreases when e increases. The spectrum of the 
synchrotron emission from these particles is also curved, with 

t^'-. (6) 

On the other hand, Ep scales like e, so a negative correlation 
between b and Ep is expected. This basically arises from a loss 
of acceleration efficiency at high energies. 

On the other hand, a connection of log-parabolic spectra 
with acceleration may be also understood in the framework 
provided by the Fokker-Planck equation (Kardashev 1962) 



dN 5/ ^dN\ \ 



dj I dy\ 



(7) 



This describes the evolution of the distribution function Niy, t) 
of the electron energies ymc^ under stochastic (first term on 
the r.h.s.) and systematic acceleration (second term on r.h.s.). 
The simple solution given by Kardashev (1962) for an initial, 
monoenergetic injection at the energy yo fnc^ reads 



N{y, t) ccy ^ exp 
with 



-(A, +A2-ln(r/ro))' 



4A, 



(8) 



Jo 



(t) dt. 



Eq. (6) again describes a log-parabolic distribution, with the 
curvature term 



1 

4aT' 



(9) 
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inversely proportional to the coefficient of the random acceler- 
ation component. 

On the other hand, the peak energy of N(y, t) is given by 

7n,a.=7oe^'-^'yoe^'-'''\ (10) 

The logarithm of the peak energy Ep of the SED is closely pro- 
portional to the position of the peak of the logarithm of Niy). 
If Niy) has a log-parabolic shape, using again Eq. (6) that ap- 
plies to any log-paraboUc energy distribution, we obtain 

3 

lnEpOc21ny„ax+Tr- (H) 
5b 

In terms of 5 ^ and Bp we see the logarithm of Ep to be inversely 
proportional to b; this is consistent with the inverse correlation 
observed between b and Ep. 

In this framework the inverse correlation constitutes the sig- 
nature of a stochastic acceleration process that broadens N(y, t) 
as to decrease its curvature while Ep is driven to higher values. 

6. Conclusions and discussion 

Using a large set of X-ray observations of Mrk 421, we have 
presented and discussed the results concerning the spectral 
variations. We have shown that correlations exist between the 
peak values Sp of its SED and the peak positions Ep, and be- 
tween the spectral curvature b and Ep. 

The former may be interpreted in the framework of syn- 
chrotron emission. The values of the power-law slopes obtained 
from our fits in Sect. 4 are bounded by a < 1 .2 ± 0. 1 . As such, 
they rule out the case a - 4 applying if the beaming factor 
5 were the dominant driver of the spectral evolution; they are 
instead compatible with a combined effect of variations of B 
(corresponding to a = 2) and of a rescaling of y (a = 1.5). 

In parallel, a secondary role for 5 is confirmed by the bound 
r < 0.3 it can contribute to the observed S p, Ep correlation 
coefficient value r » 0.7. On the other hand, we have set an 
upper limit to the beaming variance; from the analysis in sect. 
4.2 we have found a low fractional variation cr^/yUa < 10% even 
on considering conservative values 6 ^ 10. The remarkable 
implication is that the beaming factor of Mrk 421 remained 
closely constant during a time span of about 9 years. 

This limit is relevant in the framework of the internal shock 
scenario. This assumes that shells ejected from the central en- 
gine with slightly different relativistic velocities and slightly 
differing angles collide in the jet at sub-parsec scales and pro- 
duce flares. The temporal behaviour and the radiative efficiency 
of this process depend on the collision frequency and on the 
collision energetics, respectively; two versions are found in the 
literature. Guetta et al. (2004) assume that shells are ejected at 
a frequency close to 10 Hz, with F values distributed around 
the average value of about 15 after a random (uniform) distri- 
bution with a considerable dispersion, about 3. The dispersion 
is considerably larger than the values obtained from our analy- 
sis. On the other hand, Tanihata et al. (2003) assume values of 6 
following a normal distribution with erg /fig <k 0.1 and ejection 
intervals around 600 s; whence they obtain a good reproduction 
of the temporal behaviour, but also a very low radiative effi- 
ciency. The upper limit derived from our analysis, much lower 



than the value assumed by Guetta et al. (2004), emphasizes the 
efficiency problem reported by Tanihata et al. (2003). 

The correlation we have observed in Mrk 421 between b 
and Ep is interesting in the framework of the electron acceler- 
ation mechanisms. In Mrk 421 we systematically observed a 
decrease of the curvature b as the peak energy Ep increased. 
To understand this behaviour we have used in Sect. 5 a Fokker- 
Planck description of a dominant electron energy distribution. 
The solution of this equation for an initial mono-energetic in- 
jection predicts (see Eqs. 8-11) that with ongoing stochastic 
acceleration the curvature should decrease while the peak en- 
ergy moves to higher energies. A more detailed understanding 
of this dynamics requires a full computation including radia- 
tive cooUng and fixing the relative weights of the systematic 
vs. the stochastic acceleration component; this will be studied 
in a different paper (Tramacere et al. 2007 in prep.). An alter- 
native explanation of this correlation is discussed in Sect. 5 in 
terms of statistical, energy-dependent acceleration probability. 
This leads again to a correlation as observed. 

We wish to stress a conmion point to both views, that is, 
the relevance of the curvature parameter to understand the ob- 
served spectral evolution of the source. The negative correla- 
tion between b and Ep strongly indicates the dynamics of Ep to 
be related to stochastic or to statistical (and energy dependent) 
acceleration mechanisms; it would not hold with the beaming 
as the main driver of spectral variations. So the results from 
this second correlation are consistent with those from the first, 
namely, the Sp-Ep correlation. 

No significant correlation has been found between S p and 
b. This lack may arise from the opposite signs of the correla- 
tions between Sp-Ep and Ep-b adding to the considerable dis- 
persion of the data. 

Finally, another view on this matter is provided by the anal- 
ysis given in Sect. 5 (see Fig. 5) concerning the PDF of the 
spectral parameters Sp,b and Ep. It is seen that Sp and b have 
enough synmietry in their PDF to be reasonably approximated 
by a Gaussian distribution with minor deviations in the tail 
(left panels). The parameter Ep, on the contrary, shows a more 
skewed distribution (right upper panel), that could be better ap- 
proximated by a log-normal shape, i.e., by a Gaussian in the 
variable log Ep (right bottom panel). We note that a similar dis- 
tribution has been successfuUy used also to describe the statis- 
tical properties of the GRB peak energy distribution that may 
depend on broadly similar physics (loka & Nakamura 2002). 

A point to stress is that the log-normal distribution consti- 
tutes the asymptotic limit from the central limit theorem in mul- 
tiplicative form; in fact, it has been shown by loka & Nakamura 
(2002) that the limiting log-normal form is closely attained al- 
ready after 3 steps. Stochastic acceleration, for example, may 
be treated in terms of multiplication of a number of random 
fractional energy gains e, see Eq. (5). The issue will be dealt 
with in more detail in Tramacere et al. (2007 in prep.); here it 
provides complementary support to our stress on the relevance 
of stochastic acceleration to understand the spectral variations 
of Mrk 421. 
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Fig. 5. Histograms of S p, b and Ej,: filled boxes represent the XMM data set, while hatched boxes represent the FULL data set. 
The solid lines show best fits with a Gaussian, with average (ji) and standard deviation (cr) given in the labels. 



7. Appendix A: Correlations and fits in the 
presence of uncertainties on both axes and of 
hidden parameters 

We first recall the general formalism used to evaluate correla- 
tion coefficients for data with uncertainties on both axes. 



For two statistical variables defined as: 
X - {xi} - {xi + dxi) 

y = [yi) = [yi + dyi] , 



(12) 
(13) 



the covariance is expressed by 

cov(x, y) - cov(x, y)4-cov(x, dy)-i-cov(y, dy)-i-cov(dx, dy) , (14) 

and the correlation coefficient reads 
cov(x, y) 



(15) 



Here 0-^ is given by 

<Tx - cov(x, x) -1- 2cov(x, dx) + cov(dx, dx), (16) 

and similarly for ctj,. 

An issue arising in the analysis of correlations (as antic- 
ipated in §4 of the main text), is that "hidden" parameters - 
adding to the "primary" ones focused by the analysis - can in- 
troduce a covariance term, i.e., a correlation. To evaluate its 
weight we recall the general expression of the covariance term 
for functions of random variables. When we have M functions 



(/i,...,/m) of random variables xi...xn, the covariance of 
fii,fi is given by Barlow (1989) to read 

^5fk\/5fi\ 



*-«=zz(S)(||)»*--)- 



(17) 



Eq. (17) shows that even when all the variables x,- are mutually 
uncorrected (with covariance=0) a finite covariance term 
arises for the functions which share the same variables. 

We then show explicitly the formalism we actually use in 
fitting our data. We follow the approach of D' Agostini (2005), 
who discusses a development from the standard formalism (see 
Kendall & Stuart 1979) basing on Bayesian statistics. This is an 
unbiased method to perform fits on data, including uncertain- 
ties on both axes and a term of extravariance cr,, given by the 
fluctuations of a"hidden" parameter The log-likelyhood func- 
tion discussed by D' Agostini (2005) reads 



L(m,q,crv;x,y) - 



i 

(yi - mxi - 



\_ (y,- - mxj 



on having used two models, namely: 



yu, 



Xi m + q ; 



(18) 



(19) 



and a log-linear model Y = log y, X - log x with uncertainties 
evaluated on using standard propagation to yield 



Y - log Xitn + q . 



(20) 
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This method was successfully used by Guidorzi et al. (2005). 
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